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FOREWORD 


The Space and Information Systems Divisicn (S&ID) of North American 
Aviation, Inc. (NAA) under Contract AF 30(602)-3638 with the Rome Air ue 
Development Center (RADC) of the United States Air Force agreed to perform -¥ 
a 10-month study designed to develop digital computer techniques in two 
areas of interest to the RADC tracking facility. First, a differential 
correction geocentric orbit computation program for reducing observed data 
was to be prepared which would operate in a near-optimum manner at the RADC 
computer center. Secondly, a computational logic which could be utilized 
in the tracking process for driving the tracking antennae in an open-loop 
mode was to be prepared. This second program. would employ general perturba- 
tions theory in the definition of the predicted trajectory. (This former 
task is reported in SID 65-1203-1). 


This report was prepared as partial documentation of the seccnd task, 
The contents present the program logic and FORTRAN listings fcr tho main 
body of the required program. 


The program can be divided operationally into two parts; (1) the 
trajectory prediction routine and,(2) the tracking routine. 


The prediction routine uses a general perturbation theory to assess the 
first order changes in the osculating elements due to oblateness and drag. 
No singularities in inclination or eccentricity are present in the formulation & 
which is taken from the work of Anthony G. Lubowe (Appendix I). 


The tracking routine has been designed to accept as many as ten active 
stations and to output range, range-rate, azimuth and elevation data for NE 
any station viewing the satellite. 


This contract has been managed at NAA S&ID by Mr. J. A. Hill and 
directed by Mr. G. E. Townsend. J.C. Mendez, assisted by Mr. Townsend, 
designed the rationale for the program, coded the major portion of the logic, 
performed the preliminary checks of the operation, and prepared this document. 


The assistance offered by RADC personnel under the direction of 
Mr. Gordon Negus (Program Manager) is gratefully acknowledged. 
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ABSTRACT 


This document presents the formulation, computational logic, and coding 
information developed for the purpose of tracking an artificial earth satellite 
in an open-loop mode. The program was developed as a F@RTRAN IV, IBM 7094 
program which uses the standard North American Aviation monitor system (NAASYS 
version 13), The logic presented is intended for use in developing a similar 
program for the Packard Bell 250 digital computer. 


The trajectory prediction portion of the program is a general perturba- 
tion formulation developed by Anthony G, Lubowe of Bell Telephone Laboratories. 
The tracking portion of the program can accept as many as ten tracking 
stations and outputs range, range-rate, azimuth and elevation data. 
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INTRODUCTION 


This program is designed to track an artificial earth satellite in an 
open-loop mode on a Packard Bell 250 computer. Accordingly, the program 
can be considered to perform two major tasks: 


1. Prediction of the satellite's position and velocity 
vectors, and 


2. Interrogation of the active tracking stations. 
Functionally, the program consists of two driving routines: 


1. P@SVEL drives the prediction routine, and 
2. TRAK drives the tracking routine. 


Efforts have been made to insure a high level of internal consistency 
throughout the program, Particular care has been taken with the iteration 
on Kepler's equation and with the logic for computing and storing the changes 
to the orbital elements. 


The formulation for the prediction routine has been selected to eliminate 
singularities at the critical inclination, and at low (or zero) inclinations, 
and/or eccentricities; it has been taken from a paper by Anthony G. Lubowe 
which is reporduced in Appendix I. 


Briefly, the theory is as follows: a set of non-singular osculating 
elements is computed from the initial position, velocity, and time; the 
non-singular elements replace the more traditional set (a, e, i,w ,# , © ) 
which leads to the above noted singularities. The selected non-singular 
elements are: 


a=a 

=e sinw =e sin (w +a) 
v=e cosw 

p=sini sina 

q = sin i cos.n. 

FHt-@ 


BIE 


First order perturbations to these elements due to oblateness and drag have 
been included. Higher order terms, i.e., order J@, H, K have been dropped. 
Lagrange's Planetary Equation's form the basis for the oblateness perturbation 
theory. These equations express the orbit of a body experiencing a perturbing 
force in terms of the deviations of the orbital elements from those describing 
the unperturbed orbit. They are six first-order differential equations with 
time as the independent variable. The independent variable is transformed from 
time to the angle theta, in the unperturbed orbit ( © = longitude of perigee + 
true anomoly). The six differential equations can then be integrated between 
Qg and 9 by holding the orbital elements constant over the interval of integ~ 
ration. The result is a first order approximation to the changes in the elements 
due to the perturbations considered. 
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The perturbations due to drag have been taken from work by T. E, Sterne, 
A description of the work has been included in Appendix II; the description 
originally appeared in the Orbital Handbook, NASA SP-33, Part 1, Volume I. 


Input data for the program consists of geocentric position and velocity 
vectors (rectangular coordinates), the corresponding whole numbcr of days past 
zero hours, 1 January 1950, and the fractional part of a day, Due to approxi- 
mations in the formulation and the limited duration (less than 10 days) of 
most applications, no correction for motion of the vernal equinox has been 
included. Thus, the program operates in a rectangular equatorial coordinate 
system tied to the true equinox of the date of the initial conditions. Further 
required input are the step size in seconds (i.e., the interval between con- 
secutive predictions of position and velocity), the final elapsed time, the 
W/CpA and tracking station data. The main program (MAIN) reads the input data 
and drives the program. The prediction of trajectory points is accomplished 
in subroutine P@SVEL. During the first pass, the initial osculating elements 
are computed and scored; in subsequent passes, this operation is skipped. 

Then the changes to the osculating elements due to first-order oblateness and 
drag perturbations are computed. To preserve maximum accuracy, these changes 
are stored in running sums, These running sums are added to the original 
elements at each step rather than adding the changes to the elements at each 
step. The predicted position and velocity vectors are then computed using 

the updated elements, The prediction has now been accomplished and control is 
returned to the main program. 


Main next calls the tracking routine (subroutine TRAK) which computes the 
local hour angle and ‘the up, east, and north unit vectors at the tracking site, 
as well as the position and velocity vectors of the tracking site, The posi- 
tion and velocity vectore of the satellite relative to the tracking site can 
then be computed by vector subtraction, Finally, range, range-rate, azimuth 
and elevation data are calculated. At this point, control is again returned 
to the main program, 


Now, the prediction and tracking cycle is complete, and the elapsed time 
is checked to determine whether the program should continue or terminate the 
computation, Care has been exercised to assure that these operations can be 
performed in small fractions of a second on the IBM 7094 and that operation 
times on smaller computers (e.g., the Packard Bell 250) will be reasonable, 
As a result, the real-time capability required for tracking is obtainable, 


Also included in this documentation is the logic necessary to differen- 
tially correct the position and velocity vectors based on a set of observations, 
The rationale for this process invoives the reduction of a set, of observed 
minus computed residua.s; and the prediction of a corrected set of position 
and ve.ocity vectors at some future time, To be specific, the weighted least- 
squares process is used to produce agreement between a set of nine observations 
(R, Ay E, acquired at three epochs) and the corrasponding computed values at 
the specified lead time, This logic has not been included in the program 
because the acquisition and method of input pf the nine observations depend on 
the tracking and computing hardware used. It is epxected, however, that the 
observation data array 1 11d be constructed in, and the differential correc- 
tions driver routine (PRcoDIK) would be called from TRAK; thea the computed 
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corrections would be added to the computed position and velocity vectors at the 
lead time and a flag raised, When control was returned to the main program, 
the flag would indicate that a corrected position and velocity were available; 
the corresponding "corrected" osculating alements would be computed, and the 
computation would proceed as before, 


It should be noted that the lead time must be great enough to allow real- 
time computation, 


The following diagram indicates the program operation schematically, 
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SCHENATIC OF PROGRAM LOGIC 
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(PREO/CT POSITION 
ANDO VELOCITY) 
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7O POSITION 
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The following pages describe the purpose, formulation, etc. of each 
subroutine, 
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Subroutine BLZCK DATA 


Purpose: To read block data into name common 


Deck Name: BLACK 


Subroutines Called: NONE 


Functions Called: NONE 
Deck Length: 00001, 
Input/Output : 


FORTRAN | Math Common/ 
T/o | Name Name Dimension | Argument 
g | ccgn he 1 


ASTRZ 

g Ad J a ASTRO 
RE Re L ASTRG 

RP Rp ASTRG 

ALT Ay al ATMZS 

g STEP Ss 1 ATMZS 
g |DENS(M) Pm 36 ATM¢S 


ay ee 


Definition 


gravitational constant of 
the Earth 


first term in the Earth's 
gravitational potential 


EKarth's equatorial radius 
Earth's polar radius 


Lowest altitude tabulated 
in density 


distance between altitudes 
in density table 


tabulated values of density 
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MAIN ROUTINE 


Purpose: Main routine 
Deck Name: MAIN 


Subroutines Called: P@SVEL (prediction logic) 
TRAK (tracking logic) 


Functions Required: None 


Deck Length: 00515, 

Method: MAIN reads input data, calls P@SVEL and TRAK in sequence 
and exercises control over continue or terminate 
decision, 
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MAIN FLOW CHART 


Lo ig z, 
t, = t + ae 

: At 
te tr ¢ 36400 


CALL TRACK 


RETURN 
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Subroutine PQSVEL 


Purpose: This routine computes the predicted position and 
velocity vectors at time tj. 

Deck Name: POSVL 

Calling Sequence: CALL PZSVEL (RVEC, VVEC, RNEW, VNEW, 1, T1, MONK) 


Subroutines Required: ¢SCUL een osculating elements) 
computes theta for a given time) 
DBL (computes oblateness perturbation) 
DRAG (computes drag perturbation) 
VECT (computes position and velocity vectore) 


Functions Required: NONE 
Deck Length: 3128 
Input/Output: 


Dimension Description 


Comnon/ 
Argument 


current position vec- 
tor 


current velocity vec- 
tor 


predicted position 
vector 


predicted velocity 
vector 


current time 

Arg time of prediction 

Arg W= spacecraft weight 
Co= drag coefficient 

A= cross sectional are 
index = 1 for first pass 


through, 2 for subse- 
quent passes 
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Description of Equations: 


This routine predicts the satellite's position and velocity vectors at 
a specified time. P@SVEL is called each tim the position and velocity are 
to be predicted; when the prediction has been made, control is returned to 
the main program. 


During the first pass, the running sums of the perturbations to the ele- 
ments are set equal to zero. $@SCUL is called to campute the orbital elements 
and the original values of the orbital elements are stored. During subsequent 
passes, this computation is omitted and the updated elements are used to com- 
pute the current value of theta. Then, in both the first and later passes, 
the value of theta at the prediction time is computed (in subroutine THET). 
Next, the changes in the osculating elements due to oblateness (subroutine 
DOBL) and drag (subroutine DRAG) are computed and the running sums of per- 
turbations to the elements are updated. The updated running sums are added 
to the original values of the osculating elements at the prediction time. 
These values are used to compute the satellite's position and velocity vec- 
tors (in subroutine VECT). 
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SUBROUTING _PUOSVEL 


monk > / 


CALL DSCUL(R,T, ty, 2, M,Y,P)B 7) *)%) | 


Gis 
Por 
jor 
Mer 
Vor 
Tor 


n 


YX EMD O 


CALL THET (9,4, Vi %) 7, to,n) 


Cacl THET (0, M,V,9,7, tn) 


SID 65=1203-3 
Era ye 


CALL DQBL (a, PG y My Vy 7, toy € 7) 90) Os day Mp; Ig ,d ud An) 


CALL ORAG (OQ; Py Gy MV, Ada, ddp, dd, ddy dd Vv, dd, fos ty yNsie 
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a = ae tT Qsum 

P * Pr * Psur 

G = Gor * Fsum update the 
ns oe mae er orbital elernents 
V = Vor t+ Vsumy 

T= jor Tsun 


Dv 
No 
y 


F> 7T+pe 
aa 
Sum O 
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Subroutine MSCUI. 


Purpose: This routine computes the osculating elements corres- 
ponding to a given position vector, velccity vector and 
time. 

Deck Name: PSCULL 

Calling Sequence: CALL $SCUL (RVEC, RDAT, T, A, MU, NU, P, Q, TCAP, SGN, 
THETA) 


Subroutines Called: CR@SS 


Functions Called: AMAG (vector magnitude) 
DOT (dot product) 
SQRT 
ATAN2 (arc tangent) 
SIN 
cés 


Deck Length: 2h 


Input/Output: 


Definition 


Common/ 
Argument 
Arg 


position vector at time 
t 


Arg velocity vector at time 
t 


current time 


SIN SIN &% 


S/iN4 COS 
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Math 
Name 


~w 
yh 


oman) 
Dimension Argument 
ai | Arg 


J Arg 

1 Arg 

1 ASTRZ 

1 ASTRO 

aR ASTRZ 

1 ASTRO 
~ 28 « 


Definition 


F=x¢ ~n(Grn) = hime 
of equinox passage", 


= +1 for posigrade 
orbits 

= -l for retrograde 
orbits 


O=- wx 
the orbit 


= longitude in 
gravitational constant of 
Earth (length? /time2) 


first term in Jeffrey's 
gravitational potential. 


equatorial radius of the 
Earth 


polar radius of the Earth 
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SUBROUTINE @QSCUL 


/£V/— p2_q2 a an use + if OF £€9O 
B= ~ (Ff 90<2£ 4/80 


< 

z 

® 
& 


4 =Arcton (SLE 
LOS Oo 


A? ecosfsi~@ -~ EC SIN F casO 
VY = Ceosfcos@ FC SINF SING 
sie 

W-e2 © sf 


C SINE * (+ @COSF 
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SUBROUTINE @OstuL Cont) 


No 
id e> .of 


~f=- ee 2-(eswFf )* 
E-F2=-CS/ VE $cosF rE Ces 2" ] 


4 =Arctan ae) 
@ cosf 


~ 22 
SINE V7 e SIN F 


/+ @ cosf 


Cos € = cosf ee 


/¢ @cosf 


&£ = Arctan CAS) 


Cav 


7 8 7-2 Ce +E -f-eSNE) 
E 


Return 
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Description of Equations: 

Subroutine @SCUL computes the six osculating elements corresponding to 
a given position vector, velocity vector and time. 

The semi-major axis, a, is computed first. 


Let 


V =/RI 
The energy equation is 
Ag jf 2k 
Q ae r 


where k is the gravitational constant of the central body (length?/time2), 
It follows that 


PP 5. ie 
: 2k -rv2 (1) 
The second and third elements, p= :swesiw 2 and g:swi cos Ds 


are found from the angular momentum vector, B. _,The sketch shows that the 
right ascension, « , and declination, 6 , of H 


yas 
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are 


= - WF 
8y F Zz --z 


The vector, i, can be written, 


> 
H* Coos my Cos 84 1 51N Ay COS by , SINS) 


or 


~~ ‘. . . 
H:( sing SIN 2. » ~SINA COS HK, cos x) 


In terms of p and q, 


H= Cp, “G cose ) 


Now, H can be computed directly from zg and R 


> = 
7. RXR 
H = ry = CA A A,) 
/RxR | 402 ) 3 
and, finally 
ph, 


(2) 
g7 he 
The next two elements, yu:@¢s/v © ani vzecos @ , are computed 
using the variable o = & * # = "longitude in the orbit" ( @ = 
longitude of perifocus and # = true anomoly). 
Write ~ and Vv as 


M2 @snw = ASIN (O-Ff) 
= (€ cos Ft) sin @ ~(eswiticos @ 


(3) 


€ cos GW = e@ 00S (O-Ff) 


XS 
nu 


“(ecosf)cos © + (eae sinf)swe 
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The four quantities e@ecosf,eswwf 5 sw O , and cos 6 
determined below. 


(1) @ cos £ can be found by combining the equation 


re (7 -e7) 


1 +ecosf (4) 
with 
Az /Hl = kare?) 
Thus 
2 
e@cosf = ace vf (5) 
Ar 


(2) Differentiation of equation (4) yields 


_ atr-e r 
esiw f = = — (6) 
The sketch shows that * and * can be written 


r 
¥ - Vcosy 
r 


Using the expression for angular momentum 


AzrVcosy 


f becomes 


f ce (8) 
And substituting (7) and (8) into (é) 


aci-e%) (¢(R- RB) 
eswf = igs AA! 


ake 

re r 4 
Salis!) 0B ip) 2 
FA A 


- alf-eVh RR) 

rka(/~e#) 
Th ‘ 
i e@ sw f= —A— (# -#) (9) 
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And .\e first two quaniities are determined. 


As an aid in fineing s/v 6 and cos © refer to the sketch below: 


Zz 


6: orf: a+6 


From spherical trigonometry, 
A . 
SIN @ SINK = SINS 


. “A “A 
cos £ = LaNnm ~ SING 
tane cosQ sy 


N 
Cos 6 + cose cos & 


(3) Now, write syswo 


o “ 
SIN © = SIN (SutO) FeS‘nNSNL cose 7 Clos sin 6 


Substitute for a) 


SINS 


“ 
SIN @ = SIND Cosa Cosh ¥ Cos ; 
SIN 


re ae 
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A 
SIN @ = SIND. COSA COSE FCOSLLSINGA cogs § - Cos. SING COS ) 


r+cos Q Sans 


SINA 


combine first two terms and substitute for a 


SIN © = COSKH~ COS DW. cosé | Siw case] cos QL SING 


cos§ Sv dk SINE 


a 


Cos§ cosa 4#Cos asd [ Sess | 


SAN 4 


cos dcosk *+COSK SINS SINE, 
1teose 


ares (10) 
i+ COSE. 


(4) The same procedure yields 


W 


(11) 


COs OF ~ i 
r (FCOSL 


The last osculating element, t, is 
T 2 €-nCW+r) 


where M = mean anomaly and n = V a , the mean motion. 


Then, using Kepler's equation, 


F = ao 2) (arf —-F x+E& -~@ SIN E) 
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= e-n[e +(e-f)-(e@ sev &)] (12) 


The last term, esswé&, is computed from 


ee Vi~ e?(e swf) 


Jr ecost 
sre? = Vi~-py2-y? 
and @ swwF and e cos F are given by equations (5) and (9). 


If @s ¥a24¢py2 is small enough so that terms of order e* may be 
neglected &-f can be computed from the series (see Reference, page 64) 


e fe 
E-f: ~Cesmnf)fy niece cos f) + $2 aes gest) J 


If @ is so large that the series cannot be used, # can be computed 
from esswf and e cos f 


ae e sunt ) 
neeree ecosys 


Then the eccentric anomaly, E, is computed from 


Sin E =H [- €2 swf 
/+ecos f 
cos£ 2 Sut cost 


/+recosf 


And the difference, & ~¥# , can be formed. 


Reference: Brouwer, Dirk and Clemence, Gerald M., Methods of Celestial 
Mechanics, New York, Academic Press, 1961. 
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Purpose: 


Deck Name: 


Calling Sequence: 


Subroutines Called: 


Functions Called: 


Deck Length: 


Input/Output : 


Subroutine THET 
This routine computes 0 =2+G@ = longitude in the 
orbit, for a given time and set of orbital elements. 
THETT 
CALL THET (A, MU, NU, THETA, TCAP, TI, NM@T) 


NONE 


ATAN2 (arctangent ) 


477. 


FORTRAN Math Common/ 
T/o Name Name Dimension | Argument Definition 


at A 

T MU 

a} NU 

¢g THETA 

n§ TCAP 

+ 61 TL | 

g NMZT 

E GCON 


a 

Vans 
v 
6 
T 
t 


n 


wi 


1 Arg semi-major axis 
1 Arg € in te 
1 Arg @ too G) 
a Arg f+@ = longitude in the 
orbit 
1 Arg "time of equinox passage" 
1 Arg current time 
i 
1 Arg mean. motion = ; ays 
1 ASTRO gravitational constant of 


Earth = (length?/time*) 
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Development of Equations: 


Subroutine THET computes the “longitude in the orbit", 0= 4 r@& , 
given time and a particular set of orbital elements. The mean motion and 


eccentricity will be needed 


C20 urs Vv? 


If © is so small that terms of order e% can be neglected, 9 is computed 
using a series expansion. First, let & =”(¢-7F ,F=e .imM 
and @ = € coo M. 
Write § as 


Fee wnM=e sn(nt-nT) 


= € ain (Rt-aT + &-&) 


Finally 
$= UV anK — coo& 


A similar manipulation shows that 
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C= Veooks “snk 
The series expansion for in terms of M is, (Reference, page III-27) 
f=M+2e ainM + e* ain 2M + &" (/3.2in 3M-3.2inM)+... 
Substituting 
tn 2M = 2 aimMcooM 
in 3M = 3.atn M- 4 ain®M 
and combining terms 
p=M+2(e ainM) + 3 (€ 2inM)(e coo) - 


- 3(e ainM)(ecooM)*- 4 (e .aenM)° 


adding G to both sides and noting that 6 =$r @ and &= n(t-F )= M+ wo 
we have 


G=%+ Slat eros? 384] 


For the case of larger @ , 0 is computed via classic celestial 
mechanics methods. G@ and M are found immediately. 


B= Qutan(G) 


Mzep(t-T)-@ 


Then, tne eccentric anomaly, £, can be found from Kepler's equation, 


a ee SID 65-1203-3 


M=E-e¢ sine 


bv a Newton-Raphson iteration technique. £, , the first guess for E, is 
selected by truncating the series for E in terms of M 


2 
E=Mre ainM +S in (2M) 


The estimate is improved according to 


V/-e% ame 


Reference: Jensen, J., Kraft, K. D., and Townsend, G. T., "Orbital 
Mechanics, Chapter III, Orbital Flight Handbook", NASA SP~33, Volume 1, 
part 1, dated 1963. 
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Subroutine DBL 


Purpose: To compute the perturbations of the nonsingular 
osculating elements due to oblateness 
Deck Name: DOBOL 


Calling Sequence: 


DZBL (A, P, Q, MU, NU, TCAP, 


TO, TI, THETAG, THETA, 


DA, DP, DQ, DM, DN, DT, NM@T) 


Subroutines Called: PBCON 
Functions Called: SIN 
cgs 
SQRT 
Deck Length: 017508 
Input/Output : 
I/a | Name Name Dimension | Argument 
I A a 1 Arg 
I P Pp L Arg 
i Q q Eb Arg 
I MU KE 1 Arg 
{ NU Vv 1 Arg 
I TCAP t 1 Arg 
x TQ ts 1 Arg 
| E | TI | ‘, | 1 | Arg 


Arg 


Arg 


mre (ts ae 


Definition 
semi-major axis 
p =-amnl am 2. 
q =m Coo LL 
Ma= Cane 
V= eee 
"time of equinox passage" 
current time 


time of prediction 


current longitude in the 
orbit 


longitude in the orbit 
at time of prediction 
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Definition 
| change in semi-major axis 
change in p 
change in q 
change in 4 
change in V 
change in f 
mean motion 


gravitational constant of 
the earth 


first harmonic in Jeffrey's 
gravitational potential 


equatorial radius of the 
earth 
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FLOW _DB/AG RAM D@BL 


/ ay iad 92 
Vues ve 

2/RE 
a(/-e2)3 


La(/~-e2))}* 


u 
v 


CALL @BCON 


= 2dR (/t sin® +V cose) 1 -(¢ SIN O ~pcos6)* Z 
Aa ” g2(;~e2)> 4 [+ vi ] 
Op 
g 
Ap =fp {7,0 +B SINO #P,C05 8 + & sin 2o+ £5 cos 204%, sins0+ 200530} / 
{ v4 2 rs ra 
. Oo 
Ag fy 918, #8 SIN *p,COSO+ Bt sw Ze ¢ 8S cos 20+ 865/y30 + B7 C0S30 6, 
2 . 3 g g J 
0 


Ay “in {a,s + 4, SIN 8 ese ees SINZO 4 SE cos 20 + My SIN 38 


r2) 
/ / 

“M7 COS36 + 3 Uy SIN 40 = ae COS 40 + 7G Mo SIN TE ay coss0 | 

8 


t ° cat t . -_:7 re tal é. \ mrae 
aV =fn (Mer ewer )s chs $a sv20 7% YS CoSZO +MY SING 
2 


go 
4,605 30 +42y, sin 40 t2Y, COS 46 4 Yo sin seat “ucosse \' 
7 8 a 8 q7 3B 76 : 


0 
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s 3 2 
AT = 3 49 (T-4)-F7 (t-to)frr w sive, + Vcos | [! -3( SING, ~ PCOS &,) : 
2 a 


Fig) fess , [e218 +2. cin +€23 cose zr £24 SINZO +€2,C0S20 
L4+Cossk 2 2 2 


+€Ze sin30 - E27 epsse | i + Vi-e?| €3,0 + £3, SIN® ~€3,c0S 
6 6 
9 


8 
+ €345INZ0 + £35005 20+ Efe SN 39 - ae Ly ie 
i) 


t 
ae ees Le eat 2 + 6 
+ raver [E48 + VE4, SING wu E43 605 O E44S1N 2 


o 
+ E45€0S26 + E4, SIN 36 + €4,COS 36] + [ €4y sine + E4gCOS 6 
6 


+E4,,5iIN20 + E4,, COS2ZO+E4, SIN3O0+ EF,C0S 30+E4., SIN 4O 


8 

! 
+€4,¢C0ST8 + E4,sinso + ha reat 6 | | 

ct) 
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Development of Equations: 


The general perturbation expressions for the changes in the non- 
singular csculating elements due to the first order aolateness perturbation 
have been taken from a formulation due to Lubowe (Appendix I ). The 
non-singular elements 


a = semi-major axis 
p = sini sin. 
q = sin i cos i 
M = e sin Ww 
V = ecosW 
~ a 
T=t-w 
n 


have been chosen to eliminate low e and/or low i singularities in the 
perturbation expressions, and the formulation is first order in the sense 
that only terms of order J have been retained. Terms of order J2, H and 
D have been dropped. 


The development begins with Lagrange's Planetary Equations, which 
express the orbit of a body experiencing a perturbing force in terms of the 
deviations of the orbital elements from those describing the unperturbed 
orbit. These equations are six first-order differential equations with 
time as the independent variable. The independent variable is transformed 
from time to the angle theta = "true anomaly + longitude of perigee" in 
the unperturbed orbit. The six differential equations can then be 
interprated between 0, and a by holding the orbital elements constant 
over the interval of interration. The result is a first order approximation 
to the changes in the elements due to the perturbation considered. 
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Subroutine SBCEN 


Purpose: This routine computes the constants in the expressions 
for the oblateness perturbations 

Deck Name: PBCPNN 

Calling Sequence: SUBRZUTINE QBCZN (P, Q, MU, NU) 

Subroutines Called: NONE 

Functions Called; SQRT 
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Subroutine DRAG 


Purpose: To compute the changes in the nonsingular orbital 
elements due to drag over the time period ¢,-¢,.- 

Deck Name; TRAGG 

Calling Sequence: SUBROUTINE DRAG (A, P, Q, MU, NU, DDA, DDP, DDQ, DDM, 


DDN, DDT, TH, TI, AM@T, WCDA) 
Subroutines Called: DENSIT 


Functions Called: SIN 

cgs 

ATAN2 (Arctan) 

SQRT 

BESEL (Bessel Functions) 
Deck Length: 014738 
Input/Output: 


Definition 


Comnon/ 
Argument 
Arg 


Arg 


© 
ll 


semi-major axis 


= Sind SIN DQ. 
= SIN«& COS SL 
€ sin w 
V= @ cos 
change in a 
change in p 
change in q 
change in 4 
change in 
change in qT 


em OS OB BQ BQH HH HH H H 
PPP PF PP PP PP PP PB 
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Comnon/ 
Dimension Argument Definition 


current time 
prediction time 


mean motion 


WwW ,, pounds/square feet 
CoA 
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SUCBROUTINE DRAG 


a 


SIN Cl! 


“EL 
é€ 


mK 


cos w= ¥ 


G = Aretan (is as 


sin = Vp* + 92 
Cosk BWP = SIN%e 


ENTER 


+ 2 


No es 
COSL COS 82 - Peose COSKCOSSL =O 
SIN «A 
COosk sin SU = Foose Cosisin SLzO 
SIN & 
$= Aretan 608 £3/N S2 ~Y=0 
Cosnz cosfh 
Ws 95) 5 £2 
o> Ww 
2( ara 
hp = a(/-e) 
2S5/N § = = 
a a 
ee)” 
R c fe (asin & \2 
Rp J] 
hp 


Ré 
—~R 
= £00 ses lata + fp = I 


SG 
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SUBROUTINE DRAG Ccoont) 


A 


fa= DENSZIT Ap) 


Fe 1-2( 28) (re) Vase <ae 


= ae 
e:>: = 
Hp 
CD eee 
NO 
yas é je + alae A - | 
/-e2b i. 


7 (/-8e # 3e%) 
ae -o[rr ‘etetse! 
A 9 8c. (/-e2) 


Ae: 207 ai _ Brde ~ 367 

i / PRT EETT) 
AL = -OC- o [costw ee fC S a =) € 

(ie* +6e — 1S) 
rte *+ C/- et a) costwsiv & 

7 2 2 -/S ‘ 
Q $l = -DO¢/-e) {/ + gt [ae* ee, COSW SINK 
Aw = -cosk AQ 
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SUBROUTINE DRAE Ccon't) 


K 


(Ey He 
CoA | 


2n_ a*@, Fe¢ 

wy F 

(Gs 

Aa = -GC/te) tg [</-za) i,(¢) #Ze Z, C¢)] 


we 


s 76 . . 4 
Ae Zire pe @):7,:(¢) 2 § (LZ, (c) + I, (c) ] 


G = 


Mu, 


Mk = aaL (2) - 2, (c)] +c0s*w/[Z, Ce) Ze Leh am ¢ 
Of = ~K[ Zc) -2e1,Cc)] sin Ww casw 
AW = -cosk ASL 


AuzsnDhaetv(An raw) 


AV =cosMAe-L(AN+aw) 


Ap 2 COSL SINELODAL *G OD 
Aq > COSL cos DAL -~ PAD 


© 
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SUBROUTINE PRAG (cont) 


RETURN 
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Re er remrnrrmnrere aemrnme 


ET Ps 


Development of Equations: 


The expressions for the changes in the oscvlating elements due to drag 
have been based on the formulation appearing in The Orbital Flight Handbook, 
NASA SP-33 (see Reference). The appropriate pages are reproduced in 
Appendix II. In this theory, expressions are derived for the changes in 
the ciements (a, e, i, a ee), over one revolution; the following assumptions 
are made: 


(1) The density is assumed to be spherically symmetric and to change 
exponentially above perigee. 


(2) The satellite is assumed to move along the unperturbed Keplerian 
orbit for the integration range of one orbit. 


(3) The atmosphere rotates with the Earth at a uniform angular rate. 


Changes in the classic elements (a, e, 1,..,W ) can be related to changes 
in the nonsingular elements through the following relations: 


aq = Aa 

Ap = cosisinsr (AL) + q (A 2) 
aq + cosk cosa (ai) ~ p(AQ) 
AY * sin G (Ae) +V(AN +AW) 
AV= Cos ® Che) “U1 (a2+4w) 


~-(AWt A 2) 
n 


Db 
4 
tt 


These changes are "per revolution" and are multiplied by 


6, ~8o 
2% 
to yield the changes during the time period, t, ~ to: 


Reference: Townsend, G. W., ‘Perturbations,’ Chapter IV, The Orbital Flight 
Handbook, NASA SP-33, Volume 1, Part 1 (19637: 
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Subroutine DENSIT 


Purpose; This subroutine computes the atmospheric density in 
pounds per cubic foot 

Deck Name: DENST 

Calling Sequence: RHO = DENSIT (H) 


Subroutines Called: NONE 


Functions Called: AL@G (logarithm) 

EXP (exponential) 
Deck Length: 00136, 
Input/Output : 


Definition 


altitude above ret'erence 
geoid 


lowest altitude tabulated 
in density table 


distance bhctween values 
in density table 


tabulated values of den- 
sity 


DENS(M) 


pb | density at h 
- 83 - 
SID 65-1203-3 
REE LEN ET RTO TEM = an ee ee Bi Ss meee 14 : a 


Development of Equations: 


This routine interpolates a 36 point table to compute the atmospheric 
density in pounds per cubic feet. Densitiss between 500,000 feet and 
1,500,000 feet are tabulated in steps of 30,000; feet below 500,000 feet 
an error message is printed and above 1,500,000 feet the density is zeroed. 
The computation proceeds as follows: Let h be the altitude; then the 
distance, 4A, between h and the nearest lower tabulated altitude is 
computed first. The number of tabulated altitudes below h is 


m= fh - 50uv,000 4/ 
30.000 


The nunber of intervals of 30,000 feet below h is 


X= m-t 


Then 
4h=h ~ 500,000 — (30,000) x 


The density at h is computed assumming an exponential variation between 
the ™ and m#/ point 


= Om expel x22, (9 ( 4] 
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function QENSIT. 


ENTER 
es ERRCR 


4 < A, MESSAGE 


PRINT 


CALL DUMP 


m—-/] 
Ah-= h -A;. ~ (5) (X) 
RFR ere | oP joy (Smet } ] 


@ RETURN) 
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Subroutine VECT 
Purpose: To determine the position and velocity vectors corres- 
ponding to a given 9 and a given set of orbital elements. 
Deck Name: VECTT 
Calling Sequence: VECT (A, MU, NU, P, Q, THETA, RVEC, WEC, SGN) 
Subroutines Called: NONE 


Functions Called: SQRT 

ods 
Deck Length: 00402, 
Input /Ow put: 


P@RTRAN Math =| Common/ 
T/o | Name Name Dimension Argument Definition 
a A a 1 Arg 


semi-major axis 


a MU v7, i Arg M= @SIN GD 

I NU v il Are V=ecos Ww 

I P Pp L Arg P= SINK SIND 

L Q q 1 Arg Q= S/n « COSM 

i THETA 9 1 Arg |Q=/7& 

g RVEC R 3 Arg position vector 

i) VVEC | R 3 Arg velocity vector 

it 1 = + 1] for posigrade orbit 


- 1 for retrograde orbit 


gravitational congeane gf 


the Earth (Length?/Time“) 
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FORTRAN | Math Common/ 
Name Name Dimension jArgument 


Ad 


Definition 


first harmonic in Jeffrey's 
gravitational potential 


equatorial radius of the 
Earth 


volar radius of the Earth 


Development of Equations: 


This routine determines the position and velocity vectors corresponding 
to a given set of osculatine elements and 9. 


First, the position vector will be determined. The following relations. 
from the sketch and spherical trigonometry, will be needed. Z 


COS K= — 
cos 
SIN A? sin § Coss 
cos§& SIN 
siInO = 3108 
SINK 


The expression for the third component, 3 sis developed first 


A . 
= sind = SING SIV4 


setting 620 -Su 
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we Ee nod ea Oe et a ee SD ca =e. we : - ee pes 


oS SS eS 


3 


Q = sin ez [su © +e )] 
a 


SINE COS QLSINO -S/INLASINO_COS B 
and, in terms of p and q 


x? r[¢ sin 8 ~ pcos e/ 


Consider next, the first component of the position, X 


> 


+I>< 


uw 


cosx«cosd = cos (n+&)cos§ 


A A 
(cos NCOs & ~ HNO giv&) cos & 


Substituting for COSK and S/N KR 


“ 
x -[cos 1 COS8 ~ SIN. SINS C 
r C05 ¢ 


S/N QO. SINS COS: 5 
coe8 COSE sve ads 


- 


A . A 
COS 1. COS EG ~ SIN OD. COSk SING 


adding and subtracting sv sw 6, 


A A A 
- =COS2. COSO -SINIL.SING #£ SIN DLL SIN OG 


“aA 
“S/N SL Cose SIN O 


A . 
Cos(O tO) + SIND SING (/- COSL) 
SIN 4 
cos @ +$in& Siw SIN 
/¢COSE 


cose +( rreree ok 
4c0osh 


\t 
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ag te cree near retanrenrenmemeemecmmnmn a= 


(1) 


we 


next substitute equation (1) remembering that J} = swwé 
r 


x: cos O (ea G sin ~pcose | 


A = p* P \ 
X= r[( ar 40s 8 a ace ae sv] 


7AC0SL 


a similar computation shows that 


sf (1 - Sine +( h— Cos 
fo /*4COS& 


in equations (1), (2) and (3) 


r = ali-€#?) © atr-m? ty) 
1+ @cosf 1+Veose 4 sin O 


(3) 


(4) 


Equations (1), (2), (3) and (4) determine the position from a, p, q, Aw 


and 0. 


The velocity vector can be found by differentiating the position vector 


with respect to time. 


at 


“Se, 


r(2)+r6l(s - 3° )cos 8 an *3 


+#Coge 1#e0se 
2 LEt2 Sb) opera FRSC? 
o oe. \F/ Lb ad rae six / 


ee (ey ore] (Pa o 
a r(z)- [( pee (a) os di 


iad | (5) 


However, the quantities r andr must be determined. Since the position 
and velocity occur in the instantaneous osculating ellipse, the derivatives 


7 and @ are taken in the osculating ellipse. In particular 6 = 
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From the sketch 
rf * vy cosy 


where / = [R'] 


The angular momentum is 


h *|kaCr- e?) = Vika (/-pat-pe =rYcosy 
Therefore 
ro =f > Vra(l-p?-y?) (6) 
r 


To evaluate r¢ , differentiate 


a(/- e7) 


ie 
/4eCost 
to obtain 
po. rlecasf)(rf) 
eel 


a(s- ef) 


and substitute equation (6) 


—- r Cecosf) Vio Cie. yw? - 2) 
a(/-pe-y?) r 


) ecosf = f—_4 __ @ cas (0 -W) 
Aa /~ 2 - 1/2) a(/- 42-) 
finally 


Equations (6) and (7) determine 76 and fF , and all the quantities in 
equation (5) are known. 
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SUBROUTINE VECT 


QU -~? ~L*) 
[PCOS OF Kt SIN O 


7 ee ee = 
- x use tif Of KK GO 


| -/Ff Y¥0S £5 180 
/t f= p2 ge 


se = \/_he ___ (VSN -AHCOS@) 
aC/~u2~-y?) 


26= Vigad-v¥) (4) 
% = 2L(1-p8) cos +e Bsn el 
4 =2L(/-gA) siw6 t PAcos Of 
372 L4 swe ~Pcos 6] 


%22(2)-(48)L0-p8)siv@ ~g Bcos 0] 
gy 2(4) 4¢76) Lt -gA) cos 8 ~ PAsm @] 


3°32 8 (2 lg cos 6 + p sine] 


paced 


94 SID 65-1203-3 


OIJCOLIFA 
OSEOLIZA 
OFEOLIZA 
O€FOLIZA 
OCEOLIZA 
OTEOLIZA 
GOEOLIZA 
O62O0LIZA 
O82OLIZA 
OLZOLIZA 
O9ZOLIZA 
OSZOLIAA 
O92OLIBA 
O£ZOLIZA 
OZZOLIAA 
OTZOLIZA 
OOZOLIZA 
O6TOLIZA 
OSTOLIIA 
OLTOLIZA 
O9TOLIZA 
OSTOLIBA 
OXTOLIZA 
O€TOLIZA 
O2TOLIZA 
OTTOLISBA 
OOTOLIZA 
O600LI3A 
O8QOLIZA 
CLOOLIZA 
0900133A 
OSOOLIAA 
O*00LIZA 
O€OQO1L95A 
Oc00LIZA 
OTCOLIZA 


$8/81I/TG 


(iS#d 2+ 138d) © AMPAQHIAD + VBSY/(EPDSAYH#LNAY (E}DSIAA 
(iSe64 - 19825) * A9OHLE 2 WU/(2IDZAUeLOAGY (2Z}33AA 
{LI*8€d ~ AS#Ts} # LMGHLY - YAY/(TIIZAU#LOIOY (TIISAA 

($5ed - 1845) = UY {CJIIAY 
(13895 4+ ASa2d) # YY (293A 
(iSe€3 + ides} # UY (TJISAY 
¥Vaed = od 
eGe5 = €i 
Y¥ed - “T= 2d 
99nd - *T = I3 
wy / (SAND &® NOIDILYGS = LOGHLY 
(L3eNW - 1S*#NN) # CSNDD / NOT) LYOS = LECU 
INO3 / d = a9 
INO} / 0 = WV 
{2SeNW 7 1948MN 4 *T)} / JNDD = Uy 
(ViSHLISOD = 43 
(VLSHLINIS = LS 
SSOJeENIS + *T = INOJ 
(2888 - Zaad - °Y) LYS = S$so3 
(2#e#N ~- Zeek - "T) # Y = 3NO5 
du *3U “ry *NN99 /OYLSY/ ANEWKOI 
CEXIBAA *(EVIZAYU ANOISNAWIC 
SLIGYO JSIMNIOID WHOS "T- = NOS 
SLIGYO SASIMXIOITAYSLNAOD wOs °T = AIS 
*“VisHL “6 “‘d *AN *fN hy *¥Y SANSW3TS WILIEYD SHI OL DNIGNOdGS3SYUYOID 


SYDLIZA ALTIDIZA GNV NOILISOd JRL SALNdWOD ANIINOt STHL 


ON 


(NSS SISAA SS3AU SVELBHL £0 fd SAN 


- (SINSI - ANBW3LVAS Z39uNeGS NSB - 


Aw WA 


‘nw *¥) ADSA ANTLINOYENS 


LLJQZA 


#289 


SID 65~1203-3 


(3) 42 


> 


Le) oOovouUo oO 


> RPT Ree errant 


J JV 


C8EOLIZA 
OLEOLIZA 


98/81/10 


(SINGI 


ANBSW3LVLS 37uNeS 


Ndi 


QN3 
NUNLBS 


ALISA 


aaee 


SID 65-1203-3 


96 


er SS ER a 


paces 8 
ee ES 


oe 


8 NOTLISS 


u 02000 
Fs] 7 $1000 
| 21000 
Fs] 20000 
NOTLVIO1 


x — z0000 
NOTAVI01 


HLONS 1 


nda NOTLVI01 NdI Nd3 
iS; JINZONOdS3YYOI NS Nd 
$03 rl NOILI3S NIS 
G31 S3NILNOWENS 
; SINIONd AYLNS 
¥ zz000 o4 
zd ¥ £1000 13 
100% ¥ 41000 ag 
BY ¥ T1000 19 
INDI ¥ 90000 $$o2 
TOSWAS AdAL NOTIVI01 IDEWAS 
SAIGVIBVA WVYINUd GINOISNSWIONN 
3¥ ¥ 10000 rv 
TDEWAS 3dAL NOTLV5O1 TOEWAS 
10000 NIOTUO OYLSV 
os : CIWVIVVA  NOWWOI 
193A aANIANOYENS 
dV¥K FOVUOLS 


98/8T/TO 


NO 


S 


1LVrO1 


NOTLIIZS 
NOII33S 


NOTII3SS 


> ad ~~ <> > 4 


u 
q 
ddAi 


AIO VW NOWWOD 


Ad 


12000 
$T000 
€1000 
OTO0G 
s000d0 


J NOILVYIO4 


£0000 
00000 


NOITL¥VI01 


IDISAS 
LUBS 


LIZA 


Li33A 


Bene 


z 
ua 
tw 


SID 65-1203-3 


ed 
LOOHLY 
vv 

Ls 
aNd) 
WORAS 


dt 
NOI9 
TOEWAS 


97 


Soh es nearest ee td ee et dVWH 3°5VYNOLS LIDBA 
“OY 35¥G 98/8t/t0 = aes Bes eae 


SID 65-1203-3 


98 


*20500 ~ SI WLIO NI HL9N37 X9AG 


Bapragetne! nhc AOA aed tn 


Re eee 
a ord r 


ot ae 


eee haa 


Pansat. bane aed ARICA 5a - Bae 


Subroutine TRAK 


Purpose: To determine which of the several tracking stations under 
consideration is capable of observing the vehicle, and 
to compute range, range rate, azimuth and elevation of 
the vehicle relative to any visible station. 


Deck Name: TRACK 
Calling Sequence: UBR@UTINE TRAK (RDATE, VDATE, TW, TF, NUMBER) 


Subroutines Called: GHA 
UNIT 


CRABS 


Functions Called: AMAG (vector magnitude) 
DET (dot product) 
ATAN (arctangent) 


Deck Length 00501 , 
FORTRAN Math Common/ 
I/o | Name Name Dimension | Argument Description 
I RDATE R 3 Arg Vehicle position vector. 
> 

£ VDATE V 3 Arg Vehicle velocity vector. 

I TW rT 1 Arg Whole number part of Julian 
date. 

a TF T. 1 Arg Fractional part of Julian date. 

iL NUMBER N alt Arg Number of tracking stations 
considered. 

I STATN Sp 40 TRAST Tracking station data for a 


maximum of 10 stations. The 
data is arranged in vroups of 
4, 1.e., latitude, longitude, 
altitude, and station name. 


Horizon correction in case 
° . A °o ° 
horizon is not at O elevation. 


Gravitational constant of Earth. 
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FORTRAN Math Common/ 
I/O} Name Name | Dimension | Argument Description 


Ij Ad 
1 | RE Re | 
Eh RP Rp 


Description of Equations: 


First. harmonic in Earth's 
gravitational potential. 


Earth's equatorial radius. 


Earth's polar radius. 


The Greenwich hour angle and the rotation rate of the Earth are 
computed in the subroutine GHA, which is called immediately after entering 
TRAK. The following procedure is followed once for each tracking station 
under consideration. The subroutine UNIT is called to compute the 
position vector of the tracking station, as well as the up, east and 
north unit vectors at the tracking station site. Then the position of the 
satellite relative to the tracking station can be computed 


~» 


Re R -R; 


The unit relative position vector is, 


> 


a, P 
¢ Ter 


If t is the up unit vector at the tracking station site, and E is the 
elevation 


“av 


E= Aresin (u -@ ) 


The question of visibility can now be resolved, taking into consideration 
any horizon correction imposed by the reography adjacent to the tracking 
station. If E <H_, the satellite cannot be seen by the tracking station. 
Note that H., the forizon correction, will normally be zero. If the 
satellite is not visible, the computation is terminated and the next tracking 
station (if any) is considered. If the satellite is visible, the computation 
continues with the azimuth determination 


cos2 = Q-zZ 


sind 
; ‘ : : 2. ‘ 
wher. £ is the azimuth, z= is the north unit vector, and E is the east unit 
vector at the tracking station site. Then 


2 Sine 
2 Arcton aoe) 
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The Earth's spin vector is used to compute the velocity vector of the track- 
ing station. 


+ 
Sp = (9, 0, We) 


where wu, is the Earth's rotation rate computed by subroutine GHA. Then the 
tracking station velocity is 


_ a =< 


Vr > Spx Rr 


and, finally, range rate is 
> ~ + 


@ = V- Vz 


6: 1 eI 
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SUBROUTINE TRAK 


ENTER 


CALL GHA 


£2: 4 NUMBER 


SID 65~1703-3 


SUBROUTINE TRAK Ccont) 


& 


- horcer; 


cos 2 


SIND = e 
2 = Areton (S52 ) 


RETURN 
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Subroutine GHA 


Purpose: To compute the local hour angle relative to the mean 
vernal equinox, 


Deck Name: GHAN 
Calling Sequence: SUBROUTINE GHA (T, DD, GH, @MEGA) 


Subroutines Called: None 


Functions Called: None 
Deck Length: DO1NS, 
Input /Output : 


Common / 
Argument 


FORTRAN 
Name 


Dimension Description 


| 
fraction of a julian day 
in seconds 


whole number of julian 
days since 1 January 1950 


local hour angle 


rotation rate of the earth 
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where 


SUBROUTINE GHA 


ENTER 


~ _,004/ 780742 
é 1rS.2/x/0 3) d 


Gif = /00.075$4% + .9ES6473S5 d 


(2.9015 x10") fd? + WT 


Gj 


whole number of julian days past Oth 1 January 1950 
(JD 2433282,5) 


fraction of julian day 


yoo 
ad 
OQ 
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Development of Equations: 


The hour angle of the Greenwich meridian relative to the mean vernal 
equinox of epoch T is given in the Nautical Almanac as 


Amit) = 100% 07554260 + 02 985647346d 


+ (29 9015) x 10-134? + we (mod 360) 


where 
d = the integral number of days past zero hours, 1 January 1950 
t = time in seconds past zero hours of the epoch date 
«600417807417 
Ws 


a ea PR eFC P AER 
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Subroutine UNIT 


Purpose: To compute the up, east and north unit vectors at the 
tracking station, 


Deck Name: UNITT 
Calling Sequence: SUBR@UTINE UNIT (SLAT, SL@N, SALT, GHA, RP@L, RE, U, E, 
2, RT) 


Subroutines Called: None 


Functions Called: SIN 

SQRT 
Deck Lengtn: 00254, 
Input/Output: 


Common / 
Argument 


Description 


FORTRAN Math 
1/0 Name Name Dimension 


station latitude 


station longitude 


Station altitude above 
reference geoid 


local hour angle 


earth's polar radius 


earth's equatorial radius 


3 unit zenith vector at 
station 
unit east vector at 


station 


um 
U3} 
> 
Fe) 
a 


tracking station position 
vector 
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1 [cos* + Reyoe 
E 


a cos P COS (APGHA) 
= Up 1 = [COS @ SIN (APGCHA) 
U3 SIN & / 


COS (ANtGHA 
O 


(oa (A +GHA 


- SIN & SIN (A+ GHA) 


“SIN & COS (AFGHA) 
COS (A FGHA) / 


RETURN 
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Subroutine PREDIK 


Purpose: PREDIK is designed to compute a correction to the position 
and velocity vectors at a specified lead time to produce 
agreement between a set of 3 observations (acquired at 
3 epoch) and the corresponding computed values by a weighted 
least squares process, 


Deck Name: PRED 


Calling Sequence: CALL PREDIK (RAD, VEL, T, GBS, SIGMA2, CORR, RTRAK, X, Y, Z) 


Input/Output: 


FORTRAN 
Name 


Common/ 

Argument Dimension Description 

Arrays of position and ve- 
locity vectors at times 
corresponding to Ty, To» 
T3, and Ty. These data 
are expressed in the true 
| equator of date frame of 
reference and are assumed 
to be expressed in the 
units of ft and ft/sec, 


ARG An array o: times corres- 
ponding ~o the initial 
epoch (or first %, ¥), the 
three observations and the 
epoch at which the correc- 
tion to r, ¥ will be 
applied. These times are 
expressed in seconds and 


ran ha nefarnanand + anes 
New eee ww 


awe A PIOUS vo any 
| | | arbitrary epoch, 
A a 


réd set of observa 
bserved minus com- 
puted residuals). This 
vector is in the order of 
range-rate, azimuth eleva- 
tion, range-rate... etc, 
Units are ft/sec and 
radians 


oN ob 
Oo wm 
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FORTRAN Math Common / 
1/0 Name Name Argument Dimension Description 
I W ARG 


SIGMA2 W te) The weights for the 3 ob- 
servations, This vector 
must be ordered in the same 


fashion as is the vector Y, 


g CORR X(T) AnG 6 The estimate of the state 
vector at the time Te is 
obtained by a weighted 
least-squares process, 


(Units are ft and ft/sec, 


it RTRAK Pr ARG 3.553 Arrays of tracking station 
Kk position (and the corres- 
X, Y, 2 | U,E,N ponding up, east, north 


unit vectors) for the time 
of the three sets of 
observations, (To, T3, Ty) 
these data are assumed to 
be measured in feet, 


Subroutines Required: @BSERV (computes [ ay/ ax] ) 
TRANS (computes [ 8X;/ 3X.) ) 
MATMPY (matrix multiplication) 
MTINV (matrix inverse) 
Functions Required: None 


Approximate Deck 
Length: 01072 
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SYBROUT'NE PRED 


ENTER 


TIME = T(241)-T (1) 
f = RAD) 
> 

VY = YECC/) 


CALL TRANS 


RTRAK(Z) 
= RAO CIr) 
= Vez CLA) 
= x2) 

a ak ae 2 
= Z (TI) 


ym <> WD WY 


lis C9VTINV 
pene ae 
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ao el ES s = = To. SP ote re Teatro 


SUBROUTINE PREDIK  (Ctont) 


M= CO (Ts 7) 
a ee A 
= r"-rt 


CALL MIZ7INVY 
eee 


X Cis) = 9771) M7Y 


eee 


RETURN 
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Development of Equations : 


The routine will be constructed to compute a weighted, least-squares cor~ 
rection to the computed position: and velocity at a desired epoch based on a 
series of observed minus computed residuals recorded at a specified time prior 
to the epoch of the estimate. This objective will be accomplished by adopting 
a linear model which relates position errors at various epochs on the same 
trajectory (see TRANS). Thus, the errors in position and velocity at the 
times of the observations can be related to the errors at some standard epoch, 


by 


where 


tyst, = the epochs of observation and reference, respectively. 


Further, the errors in the observations can be expressed as linear func- 
tions of the error vector (X) (assuming that X never becomes large). This 
step is accomplished as follows: 


Y; = M vector of observed/computed residuals 
= By Xy 
where 
B= o% = Mx 6 matrix 
t ox; 
Thus, 
le Sead 


i? By (tis to) x 


Now the total set of observations can be expressed as 


Y 8, 
—_ B -_ a 
V2 Ye 2 2% ee ee Ae 
Yn 8, 
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ams enn ETD re een eee anIEnnamansinneteeanst nena la SERIE R otal anics doen Pate. Fach es oa Sn, eres Peers 


where 


“<f 
ul 


mn-vector 


oO 
as 


mn by 6 matrix 


Finally, the errors at the time at which the correction is desired can be 
computed as 


¥(T)= p(T, 4,) ¥Cto) 
or 

—_—_ =f ~<A 

Y=CP(7t,) X(7) 

2M X (7) 
The problem thus reduces itself to the derivation of a computational alogrithm 
which will construct an "optimum estimate" of X(T) for the case where mn >6. 
Assume that the vector Y is now contaminated with noise as: 


= 


Y 


MX(T) + 
or 


2 = Y-MX (7) 


m 
Further, construct the comparison function F (F = w,; © ) which is desired 
to be as small as possible. al f 
Fo 1 2 
where 
W = diagonal matrix at weights 
] 
2 0 
ON 
ist aN 
~ 
ed 
O ane 
Cen 


Now differentiate F with respect to xX 


AF? ~A% (7) M™w(9~ Mk (7) -(- MET) WMAY (7) 


=-[(m"’ wy - Mw X (7) AX (t)] '-[(M" wY-M"wWM Y(7)) ax (2)] 
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Thus, if F is to be a minimum (1.e,, SF = 0), the function 
Mw ¥ = MT wom X(T) = 0 

and X(T) may be obtained as: 
Rt) = (Twa)? MT wy 


It is important to note that since M is dimensioned mn by 6, neither it nor 
its transpose is invertible. Thus, the equation cannot be further simplified, 


The actual problem will be a specific case of this analysis in that only 
a slightly over-determined set of equations will be processed, That is, only 
a slight amount of smoothing will be employed. In this case it is assumed 
that the following data are available, 


Tracking Station 


Up, East, North 


where : 
PR range rate 
Y.z¢ AA azimuth 
AE elevation 


2 
f Tap variance in range-rate for ith observation 
2 


variance in azimuth for ith observation 


2 2 o cJ 
k: Che variance in elevation for its observation 
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Subroutine @BSERV 


Purpose: QBSERV computes the 3 by 6 matrix of partial derivatives 
of the observables with respect to the state for the 
case where range-rate, azimuth and elevation are 


acquired, 
Deck Name: BSN 
Calling Sequence: CALL GBSERV (RVEC, VVEC, @BSN, RTRAK, X, Y, Z) 


Input /Out put: 


FQ@RTRAN 


Common/ 
Argument 


Dimension Definivion 


The positiun and veloc 
ity vectors for the 
Satellite on the esti- 
mated trajectory in the 
true equator of date 
frame of reference, 
Units are ft and 
ft/sec. 


<I 
w 


ARG 


5 
VVEC 


g BSN B 3 x 6 ARG The matrix of partials 
of range-rate, azimuth 
and elevation with 
respect to errors in 
position and velocity 
(state), Units are 
ft/sec per (ft, ft/sec) 
and rad per (ft, 


f+/sec), 

Tha aaneti fan eneenn Le 
ae Ppwos 240i VEULUL av 
the tra <ing station at 


the time of the obser- 
vation being considered 
(ft), 


The up, east, north 
unit vectors at the 
tracking station, 
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Subroutines Required: CR@SS (crossproduct ) 
Functions Required: AMAG (magnitude of a vector) 
SQRT (square root) 
D@T (dot product) 


Approximate Deck 
Length: 003355 
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-» ee ee a ae a SRP es ~tor per mememner mn ewer = 


we 


SUBROU T/NE OBSERY 


Partials of ronge-rate 


» 


aA 
=P ' Vp 
7 — 
= Me — 
wc’ 


£e 
lad 


Partials of AzsmuthA 
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SUBROUTINE C&SERY (cont) 


Partials of Elevation 


RETURN 
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Formulation: 


@BSERV is constructed to define the partial derivatives of the observables 
(i.sec, range-rate, azimuth and elevation) with respect to the state vector 
(i.e., ar, 4¥), This information is to be presented in the form of the 
3 by 6 matrix illustrated in the following equation: 


AR 
AA >= BC) &xce) 
AEX 
The matrix b(t) will be constructed from the following equations 
. A w™ 
ei a vy) 
A = tan” ( ren 
FA 
= -¢ —_ “ 
EL * sin (Bae) 
Cash 
pe ae 
™ _ 
ns = t | R 
_~ —_ 
Vp Vv - v, 
where: 
) 
R = range rate 
Az = azimuth 
EL = elevation 
NA 3 . : 
U,E,f = up, east, north unit vectors at tracking station 
— 
r = vehicle's position in equatorial frame of date 
rh = nominal position on reference orbit 
P = tracking station's position vector 


when it is noted that for the purpose of differentiation, the nominal position 
vector and the tracking station position vector are constant, i.e., 


This set of operations has been performed, and 
are presented below: 


the results of the analysis 
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PEN E> : a a ae ene. ese Awketeretneg, Nhe - wew-G. a TT AS ae Ce PET chi OIE CLE EE ee 
~ ed a 


l. Partials of range-rate 


Pam Abe ary Zzrs 
BOR et a a ay 
2k .( 98 aR) _ Xp RX Yr _ AX Zp Rz X y, zZ 
a CRG) GE BE By Fe Bn BB) 


2. Partials of azimuth 
= 2 ~ 2 
572 (FE) + (ER) 


3. Partials of elevation 


— A =n 
QEL | otf oat ) (Gee), ey (Ee), 
ay ‘lar ae / "| ns te 
Us ~ Z,. a7 O,°;0 
S , R*s z 
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Pcoreted 


INPUT DATA 


All input data is read in under a floating point format. The first 
card contains the initial position and velocity vectors (see Figure 1), 
The first two fields of the second card refers to the initial time. The 
first field contains the whole number of julian days past Oh 1 January 1950 
(JD 2433282.5), and the second field contains the fractional part of the 
day. The third and fourth fields contain the step size and the final 
elapsed time (in seconds). The fifth field contains the W/CpA in pounds per 
square foot, and the last field contains the number of tracking stations 
considered (input as a floating point number). The third card contains data 
describing the first tracking station. The first half of the first field 
(first six columns) holds the station name (6 letters maximum), and the 
second half is blank. The next four fields contain the longitude, latitude, 
altitude, and horizon correction for the sta.ion. The last field is blank. 
A maximum of ten tracking stations may be active; data for each active 
station must be input on a separate card following the first station data 
card. 


139 SID 65-1203-3 


UCTYBYS YOR IOJ pres auc 


Byep UuoT}BYS ZuTHoes] 


(Sep) suoTyOSII0D UOZTIOY 


€Sbe IeAo apnataTer 
! 
sap SPNyTIETE. | 


= 
= 
= 
a. 
= 


(OT = “xeu) suotzeqs 


im ~ — yooy erenbs [ spunod ut wd9/mi 
——$s <$ _—__ —c_ - ———_ 
2 Spuooes UT out, pesde[Ts [euTy 


‘Spuoo0es UT ozTs 


UTIOeLIY FO Toqumu| 


Xe Sl9zeT 9) UOTIeYS FuTyoery Fo sureuf. 


Oo Vi ¥ @ 


€ 


nw 


3 08 eee were wees vewe creer overs 


aye 


kep © jo yred ~TeuotzoRss aut} 


OS6T ‘UBT yO Wed Shep Jo JequmU eToum| Tet tut} 


Z 
Sake ee 5 
ZOJPOSA AYTOOTOA TeTYTUT JO qusttoduos x 


Z 


O08 


ZOJOOA UOT}ISOd [eT ppuy FO JUSUOMUIOS xf 


HONMd AJY LON OG NOlidIYDS3AG 


: ~"ON BOr 40 39Vd 31V0 


ViVd WWID3S0 LIDIG Ol 


NOILVOISILN3G! 


G3X1!4 


YSWWVYD08d 


NVYLYOS 


(TTPAY Were APH TTT 


YAESWNN 


~"ON W230 


SID 65-1203-3 


140 


OUTPUT DATA 


Immediately after reading the input data, a record of the input data is 
printed out for reference, The initial osculating elements are also printed 
out. 


Beyond this, output consists of the position and velocity vectors at each 
step and tracking data at any step when the satellite is viewed by a tracking 
station, 


To illustrate the output format and te assist in the checkout of any pro-~ 
gram developed from this document, several pages of sample output have been 
included, The initial conditions nave been printed out; note that the step 
size was four minutes, 
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SAMPLE PROBLEM 


A sample problem was constructed around ECHO II tracking data recorddd 
at Floyd Satellite Communication Terminal, Floyd is located at the Rome 
Air Development Center in Rome, New York; more precisely, the station's 
coordinates are; 


N= 284.6596 (aast) 
@= 43.1972 (north) 
h= 164 km 


The osculating orbital elements for ECHO II at epoch O min (U. T.), 
0 hours, 20 April 1965 were provided, 


a = 7528.31 km 
e= O24h7 
i= 81,45° 
w= 34.156° 
= 31,2020 
M= 113.0929 


and these were used to compute the position and velocity vectors at 
epoch in the geocentric coordinate system (true equinox of date). 


59.5, 94,38 
2 =|-2918. 1582) kan 
3783. O742 


4 (72-Th4h 4520 
V =-2.729 9969] km/sec. 
-6.07h 8353 


The position and velocity vectors were input into the general perturbations 
program ( refer to Figure 2 ) and subsequent positions and velocities 
were predicted for a period of seven days at steps of 6480 seconds, 
approximately one revolution. During the eighth day, 27 April 1965, 
ECHO IT was sighted several times by Flovd Tracking Station in New York. 
Consequently, at the beginning of the eighth day the step size in the 
program was cut to 60 seconds. After each step, the program tested 

to see if ECHO II was visible from Floyd Tracking Station. If the 
satellite is seen, the program outpute range, range-rate, azimuth, and 
elevation data. This computed data was compared with the observed data 
provided by RADC; the results for passes 6069 and 6070 are plotted in 
Figures 3 and 4, 
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Results of Sample Problem: 


The Figures 3 and 4 present plots of observed and computed azimuth and 
elevation for the Echo II pass occurring at approximately the 16th hour of 
27 April 1965, The curves show that during this eights day of prediction, the 
computed position logs the observed position by 20 to 30 seconds, Using a 
speed of about 8 km/sec, a corresponding position error of 160 to 240 km is 
indicated, The angular difference between the observed and predicted azimuth 
and elevations is on the order of ten degrees, 


The time discrepancy of about 30 seconds after some 100 revolutions is 
approximately one-third of a second per revolution, The orbital period of 
Echo II is approximately 108 minutes (6480 seconds) so that the time discrep- 
ancy per period is one in the fifth significant figure. Since the initial 
osculating elements were to four or five significant figures, a difference of 
the above magnitude is to be expected, 
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ACCURACY TESTS 


A check run was made against an Encke-type integrated trajectory generated 
by a North American, SID, digital computer program (AP-110)., The comparison 
was made over three revolutions with an interval of four minutes between pre- 
dicted points, To test more accurately the real-world effectiveness of the 
prediction program, all of the perturbations included in AP-110 were activated, 
Specifically, the oblateness through the fourth-order harmonics, drag, lunar 
and solar gravitational potentials were included, Because AP-~120 predicts 
position in nautical miles, comparisons are mada in these units, The table 
shows the position results after one, two and three revolutions, (G.P. denotes 
the general perturbation formulation. ) 


initial aonditions: 


AP-110 
X 3600,0000 3600,0000 
Y 0.0000 0,0000 


Z 0.0000 0,0000 


after one revolution: 


X 3596,53910 3596,53840 
Y “124 ,31783 -124, 25333 


vA 102,55489 -102.47667 


after two revolutions: 


x 3586, 24040 3586,17630 


b § -248,01990 ~248,14364 


Z -204 60541 -204,69866 
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after three revolutions: 


x 3569,17900 3568 .,97940 
Y -370,89563 -371,66846 


vA -306 05043 -306.77181 


The results indicate agreement with the integrated "real-world" trajectory to 
within one nautical mile after three revolutions. 
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aa CONCLUSIONS AND RECOMMENDATIONS 


The accuracy test with an integrated Encke-type trajectory and the results 
of the sample problem have confirmed the soundness vf this approach to the 
task of driving a tracking antenna in an open-loop mode, 


4 The scope of the stody and the time restrictions have precluded an exten- 
4 sive checkout that would elevate the status of the F@RTRAN IV program to the 
“a all-inclusive level of "operation," Further checks on the variation of 


accuracy with prediction step size along with a complete incorporation and 
4 checkout of the differential corrections section would be necessary before the 
. status promotion could be made, 


y Specifically, the following items merit consideration: 
1, Increase program efficiency by combining groups of subroutines, 


4 2. Increase program flexibility by including the effects of solar 
04 radiation pressure and lunar and solar gravitational potential. 


‘ 3. Provide for initial-condition options, e.g. the non-singular 
osculating elements and/or a more conventional set of orbital 
elements, 
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